# Solving systems of equations
# To solve a system of equations, we will use the Python module
# 'sympy' and the function 'linsolve'. We first have to tell Python which
# variables we will be using. In this first trivial example, we will use "x"
# as the variable.
# Import sympy:
from sympy import *
# Define our variable:
x = Symbol('x')
# Next, we define our equation(s). It is assumed that the expression is
# equal to zero. For example, the RHS of the equation is x - 2 and the LHS
# is zero.
eqn1 = x - 2
# Now we can implement linsolve([...], [...]). In the first set of
# square brackets, list the equations. In the second set of square brackets,
# list the variables that you want to solve for.
soln = linsolve([eqn1], [x])
soln
# Python returns the solution as a `FiniteSet'. We can covert the FiniteSet
# to list via:
solnList = list(soln)[0]
solnList
# Finally, we can index the list to get the final solution
solnFinal = solnList[0]
solnFinal
# Here's a less trivial system of two equations and two unknowns.
# First, we need to define a second variable.
y = Symbol('y')
# Here's the system of two equations and two unknowns.
eqn1 = 10*x - 3*y - 5
eqn2 = -2*x - 4*y - 7
# We now implement linsolve.
soln = list(linsolve([eqn1, eqn2], [x, y]))[0];
print('x =', soln[0], 'and y =', soln[1])
x = -1/46 and y = -40/23
# It is also possible to algebraically solve a system of equations in terms
# of other variables. Below is a system of 6 complex equations and 6 unknowns.
# The example below solves for the currents in an all-pass filter.
# (Optional experiment in PHYS 231.)
# We must define the six currents that we wish to solve for as variables as
# well as the other relevant variables (R is resistance, C is capacitance,
# L is inductance, v is voltage amplitude, and w is angular freqeuncy.
v = Symbol('v')
R = Symbol('R')
w = Symbol('w')
L = Symbol('L')
C = Symbol('C')
# Here is a method to define all of six of the i0 to i5 symbols in a single line
i0, i1, i2, i3, i4, i5 = symbols('i0:6')
# Before we write our system of equations, first note that in Python you make
# a number complex by following it with a 'j' (no space). For example,
# to enter z = x + j*y we would type 'z = x + y*1j'. In these expressions
# j represents sqrt(-1). Below we evaluate the square of of j which should
# be equal to -1.
1j**2
(-1+0j)
# To continue an input on a new line, use the backslash notation.
eqns = [i0 - i1 - i2, i1 - i3 - i4, i2 + i4 - i5,\
i0*R + 1j*w*i1*L + i4*R + 1j*w*i5*L - v,\
i0*R + i2/(1j*w*C) + 1j*w*i5*L - v, i3/(1j*w*C) - i4*R - 1j*w*i5*L]
# For some reason, 'linsolve' doesn't work here. However, we can use 'solve'.
soln = solve(eqns, [i0, i1, i2, i3, i4, i5])
soln
{i0: (I*C*L*v*w**2 + 2.0*C*R*v*w - I*v)/(2.0*I*C*L*R*w**2 + 2.0*C*R**2*w + 2.0*L*w - 2.0*I*R), i1: I*v/(-2.0*L*w + 2.0*I*R), i2: -C*v*w/(-2.0*C*R*w + 2.0*I), i3: C*v*w/(2.0*C*R*w - 2.0*I), i4: (-I*C*L*v*w**2 - I*v)/(2.0*I*C*L*R*w**2 + 2.0*C*R**2*w + 2.0*L*w - 2.0*I*R), i5: -v/(-2.0*I*L*w - 2.0*R)}
# As you can see the solutions is very complicated!
# The output of 'solve' is an object called a 'Python dictionary'.
# The individual solutions can be accessed using soln[i0], soln[i1], ...
# Here is the solution for i0:
soln[i0]
# We can use 'subs([],[],[],...)' to enter in numerical values for the various
# symbols. I'll import 'numPy' too so that we can access pi. Here is a numerical
# value of i0.
import numpy as np
i0num = soln[i0].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)])
i0num
# We can force Python to tidy up the number using 'N()'.
N(i0num)
# We can use an option with N() to specify how many sig figs to keep. In the
# output, the numerical value shows only 1 sig. fig. because it is exact.
N(i0num, 3)
# Here are numerical values for all six of the currents.
print('i0 =', N(soln[i0].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i1 =', N(soln[i1].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i2 =', N(soln[i2].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i3 =', N(soln[i3].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i4 =', N(soln[i4].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i5 =', N(soln[i5].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
i0 = 0.000500 i1 = 0.000194 - 0.000244*I i2 = 0.000306 + 0.000244*I i3 = 0.000306 + 0.000244*I i4 = -0.000112 - 0.000487*I i5 = 0.000194 - 0.000244*I